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Abstract 


In this paper, some monotonicity-preserving (MP) and positivity-preserving 
(PP) splitting methods for solving the balance laws of the reaction and 
diffusion source terms are investigated. To capture the solution with high- 
accuracy and resolution, the original equation with reaction source term 
is separated through the splitting method into two sub-problems including 
the homogeneous conservation law and a simple ordinary differential equa- 
tion (ODE). The resulting splitting methods preserve monotonicity and 
positivity property for a normal CFL condition. A trenchant numerical 
analysis made it clear that the computing time of the proposed methods 
decreases when the so-called MP process for the homogeneous conserva- 
tion law is imposed. Moreover, the proposed methods are successful in 
recapturing the solution of the problem with high-resolution in the case 
of both smooth and non-smooth initial profiles. To show the efficiency of 
proposed methods and to verify the order of convergence and capability of 
these methods, several numerical experiments are performed through some 
prototype examples. 


AMS subject classifications (2020): 65N08, 35L65, 65120. 


Keywords: Balance laws; Splitting method; Monotonicity-preserving. 


*Corresponding author 

Received 26 August 2020; revised 22 October 2020; accepted 12 November 2020 
Fayyaz Khodadosti 

Department of Applied Mathematics, Faculty of Basic Sciences, Sahand University of 
Technology, Tabriz, Iran, e-mail: fayyaz64drQgmail.com 


Javad Farzi 
Department of Applied Mathematics, Faculty of Basic Sciences, Sahand University of 
Technology, Tabriz, Iran, e-mail: farzi@sut.ac.ir 


Mohammad Mehdizadeh Khalsaraei 
Faculty of Mathematical Science, University of Maragheh, Maragheh, Iran, e-mail: 
Muhammad.mehdizadehQ@gmail.com 


73 


74 Khodadosti, Farzi and Khalsaraei 


1 Introduction 


Balance laws appear in many different areas such as fluid dynamics, gas 
dynamics, chemical reactions [13, 19, 1, 14, 29, 35, 37]. These equations are 
hyperbolic and are sometimes called inhomogeneous conservation laws in the 
literature on physics. In particular, these equations are described as follows: 


(1) 


uo(x) = u(x,0), « ET, 


where g(u) and f(u) are smooth functions of u and J is the computational 
domain. The characterization and design of high-order and high-resolution 
schemes is a challenging task in the numerical analysis [12, 19, 1, 24, 25, 26, 
29, 4]. The main goal of this paper is, therefore, to design schemes with 
high-accuracy and resolution to obtain the true solution of the balance laws 
including reaction and diffusion terms [19, 14, 7]. One way to deal with such 
problems is to split the original problem into two sub-problems and then solve 
each problem independently by applying an appropriate technique, which is 
sometimes called the fractional step approach [5, 28, 29, 32]. To this end, the 
problem (1) splits into the following sub-problems 


P, 2 Ut + g(u)« _ 0, 


Po: uw = f(u). 


In other words, the original problem (1) splits into the homogeneous conser- 
vation law P; and the ODE problem P:. Mathematically, each sub-problem 
is computationally solved by a standard numerical scheme. However, solving 
the initial value sub-problem P; by a suitable numerical approach such as 
high-resolution shock-capturing methods, and then considering its solution 
as an initial profile for the sub-problem P, will eventually lead to the solu- 
tion of the original problem as we solve sub-problem P2 with an ODE solver 
such as Runge-Kutta (RK) or multi-step methods. For a semi-discrete dis- 
cretization such as linear multi-step methods, we can study the global error 
propagation [10]. This procedure is regarded as a simple method to deal 
with this problem and often produces good results. Also, using this process 
the high-resolution methods such as the total variation diminishing (TVD), 
weighted-ENO (WENO) and MP methods can be applied directly for solving 
the homogeneous conservation law P, [1, 29]. 

In the past few decades, some numerical methods for evaluating the solu- 
tion of the balance laws have been applied. For instance, Strang [4] proposed 
a splitting method for solving the balance law equations. Hundsdorfer and 
Verwer in [19] investigated one-step and multi-step methods using a splitting 
method to solve the advection-reaction equations. Moreover, they analyzed 
and extended the operator splitting method for a wide range of balance laws 
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including advection-diffusion-reaction equations in one- and two-dimensional 
spaces [1]. LeVeque in [29] introduced an analysis of numerical methods for 
solving balance equations. In [23], the corrected operator splitting method 
for solving nonlinear parabolic equations of a convection-diffusion type has 
been presented. Langscth et al. in [28] studied the order of convergence 
for operator splitting for balance laws. In [35], an A-stability achievement 
for operator splitting type time integration methods applied to advection- 
diffusion-reaction equations with possibly indefinite source terms has been 
achieved. Nuri et al. obtained the numerical behavior of a modified Burger’s 
equation using cubic B-spline collocation finite element method after split- 
ting the equation with Strang splitting technique in [40]. Holly and Polatera 
in [18] presented a numerical method to model the contaminant dispersion in 
two-dimensional tidal currents. In fact, they applied the high-order bi-cubic 
Hermite interpolation with characteristics to obtain the true solution in the 
presence of the advection part as well as the well-known Crank-Nicolson 
scheme for the diffusion part. One of the methods of discretization of the 
homogeneous conservation laws is the use of MP methods. Some numeri- 
cal schemes with the monotonicity-preserving property are provided in the 
references [15, 5, 9, 2, 15, 17, 31, 17]. Suresh and Huynh [5] proposed an 
accurate MP scheme to homogeneous conversation laws. In the MP scheme, 
starting with an original numerical flux computed by any high-order scheme, 
a local interval is computed by enlarging the first-order MP interval. Then, 
the original numerical flux is preserved or replaced based on the relation 
between the original numerical flux and this local [8, 6, 16, 2, 13, 5]. This 
process is, however, capable of preserving the monotonicity property and, in 
contrast to conventional TVD methods, maintains the accuracy of solution 
in the extremum points. The MP scheme has a simple algorithm to achieve 
a high-order of accuracy. Moreover, this scheme is widely applied to solve 
hyperbolic partial differential equations [1, 6, 2, 15, 16, 5]. In this paper, 
we focus on designing the MP splitting schemes based on the idea of the 
MP scheme for solving balance laws. This process has some advantages in- 
cluding preserving monotonicity and positivity with high-accuracy, having a 
high speed of implementation due to the MP process [5], higher resolution at 
discontinuous points and yielding accurate results in the smooth region. 


The outline of this paper is as follows. Details on the construction of 
the operator splitting schemes are presented in Section 2. In Section 3, the 
construction of the MP process for homogeneous conservation laws will be 
described in detail. The main results for the MP and PP splitting process 
for balance laws will be discussed in Section 4. To examine the constructed 
methods, some numerical experiments are given in Section 5, and finally, 
some concluding remarks are given in Section 6. 
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2 Operator splitting schemes 


Here, we discuss the discretization of the problem (1) corresponding to sub- 
problems P; and P2. For this reason, we consider two types of the operator 
splitting schemes. The first type is the so-called basic splitting scheme and 
is applied to the sub-problems P; and P2 as follows: 


d 

ae te) = L(t,u*)(t), for ty <t<tn4i, u*(tn) = Un, (2) 
d 2K 2K OK ae 

ae Q=F(0")@) for te <i Sta, ee) aw Ga), (3) 


where u, stands for the approximating u(tn), L(t,u(t)) denotes the dis- 
cretized homogeneous operator for conservation laws and F’ (t, u(t) is the 
operator for reaction source term. Finally, the solution of the problem (1) is 
obtained as we set u”*+ := u**(tn41) in the interval [tn,tn+1]. Also, the basic 
splitting scheme for the sub-problems P, and P, provides an error bound for 
the solution of order O(At) [19, 4], where At is the splitting time step. The 


second type is the Strang splitting method and is implemented as follows: 


d 


qe SL Oy tye tie 


uw (tn) = Un; 


d 

qe © =F(t,u")@), tn<t< tro, w* (te) = U* (trys) 

P (5) 
qe =L(tu™")\(), tae 6X thai, (tay) =U (Enea); 


(6) 
where the operators L(-,u(-)) and F(-,u(-)) are the same as those in basic 
splitting scheme. The final solution of the problem (1) is obtained by setting 
ut! = u**(t,11). Strang splitting method introduces an order O(At?) 
[19, 29, 40]. The homogeneous conservation laws in the sub-problem P, 
are discretized via step-size Ax = Tip —%j;_1 na uniform grid and is 
defined by the points x; = jAx, j = 1,2,...,.N, with cell boundaries given by 


,.1 = 2; + 4” in the conservation form 
Jj+3 J 2 


du;(t 1 
Be pt 5)) 


where fj+ 1 are numerical fluxes. Our aim is to calculate ce at time level 
tn4i- It is worth noting that the algorithm for solving equation involving 
time is based on the well-known three-stage, third-order TVD-RK3 method 


[15, 36, 5]. Consider : jt4 as the original numerical fluxes which are approxi- 


mations of the point value u(x 54 1, t) at the time level ¢t,,. We use an original 


numerical flux with a five-point stencil. Since, we need at least five points 
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Table 1: The constant coefficients 8; for k = 3,4, 5. 
k B5 Ba B_3 B 2 Ba Bo By Bo Bs Ba Bs 


3 1 5 101 319 107 19 1 


140 84 420 420 210 210 105 
4 1 Al 199 641 1879 275 61 11 1 
630 2520 2520 2520 2520 504 504 504 504 
5 1 61 703 371 7303 20417 15797 4003 947 17 1 
2772 13860 27720 3960 27720 27720 27720 27720 27720 3080 2310 


to distinguish between extremum and discontinuity points. Therefore, one 
choice is to consider the original numerical flux of order O(Az*) as follows: 


: 1 
faa = 5p (291-2 = ilgg ate PO Tea = 3gi42)) (7) 


Furthermore, to increase spatial accuracy, the original numerical fluxes of the 
(2k + 1)-th order in the following form 


k 
Fi+n = S- Biggj4ts k= 3,4, 5, 
l=—k 


can be used. The constant coefficients $8; are shown in Table 1. Moreover, 
compact numerical schemes are relatively high-order with a small stencil. 
However, these numerical schemes can be applied as another choice [15]. This 
process produces some oscillations in the numerical solution as the above- 
mentioned numerical fluxes are assembled. To overcome this difficulty, we 
will discuss the monotonicity property for homogeneous conservation laws in 
the next section to obtain accurate numerical schemes inspired by the process 
of MP as proposed in [5]. 


3 Construction of the MP process for conservation laws 


In this section, we will study the construction of the MP process for conser- 
vation laws with numerical flux (7). To this end, the numerical flux (7) is 
considered to be the original numerical flux Faz: Because of this, some oscil- 
lations near the discontinuous points are produced as we utilize the numerical 
flux ce 41. Therefore, we limit the numerical flux endowed with MP process 
in such a way that the oscillations near the discontinuous points are damped. 
To construct an MP process for the numerical flux (7), some definitions are 
required. The first definition is the minmod function for n arguments as 
follows: 


minmod (21, Z2, tie) = s. min ( | 21 |,| 22 |, | 2n | y 
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where 


s= 5 (sigm(21)-+sign(2s) x 


5 (sigm(21) + sign(2s)) a (sign(a) + sign(n)) | . 


The second one is the interval [ min(21, 225 +++) Zn), MAX(Z1, 22, ---, Zn)| denoted 
by Iz, ZQ5 0005 Zn] - The third one is the local curvature 


MM 
diya = minmod (dj, d;+1), 


where 
d; = Uj4+1 — 2u,; + Uj—1- 


For simplicity of notation, we omitted the superscript n, e.g., uj denoted by 
u,;. We consider the original numerical flux of five-order in the formula (7). 
Then, this original numerical flux is maintained or replaced according to the 
following limiting procedures. More precisely, if we suppose u; < f; 44 < 
uj+1 some simple calculations will show that the derived solution, ame lies 
somewhere between u;—, and u;, which yields MP-property. Therefore, these 
ensure that f;,1 lies between uj; and u¥” as follows: 
2 
ul™ = uy + o(uy — uj-1), 

where a > 2. According to the above-mentioned assumptions, we derive a 
first-order MP interval I[u;,u“?], where wu”? is defined in the following 
form 


MP _,. . ee ae 
u = uj +minmod(uj41 uj, a(u; uj-1)): 


For simulation, the condition for which the original numerical flux f +3 


uMP) 


lies in I[u;, , is equivalent to 


(hag — 1) (rg 0) <a, 


with a tolerance p = 107'°. For Fin € Ifuj,u™?] the accuracy of the 
MP method decreases to the first order in the vicinity of the extremum 
points. To increase accuracy in the vicinity of the extremum points, the 
intervals I[uj,uj41] and I[uj,u/"] were enlarged as I[uj;,uj4i,u/?] and 
TlmgeuP aa respectively, where 


, 1 aay 
uM? = 5 (uj + sti) — 544% 
and 4 
ul? — ~(3u; —uj_1) + <a 
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Two intervals T[uj, ule, | and T[uy, Uj415 | enlarge only for the 
local non-monotonicity numerical data, and degenerate to I [uj, uw? "| and 
I{u;,uj41] for the local monotonicity numerical data [5]. In practice, we 
. MM M4 
replace d; +4 with d, i and define 

M4 . 

dyya = minmod (4d; i dyj41, Adj+1 — d;, d;, dj41) . 

Intersection I[u™™, u™*] of the two intervals I[u;,uw7",u”°] and I[u;, uj41, uv? 
can be computed by 


um — max [ min (uj, Uj+1; uM?) min (uj, ue ul)! F 


u™e* = min | max (uj, Uj+1; uM?) max (u;, he ul) ; 


which indicate the local interval I[u™™, w™?*] for the MP process. The nu- 
merical original flux is replaced by the nearest bound of the interval if the 
original numerical flux is outside of the interval. Therefore, we have 


fiaal = ies +minmod(u™*— f,,.,ui" — fi.) 


Summarizing the above scenario, accuracy and monotonicity are preserved, 
which means the MP process is of high-order and high-resolution in the case of 
the smooth and non-smooth initial profiles, respectively. In comparison with 
the other existing methods such as the fifth-order weighted-ENO (WENO5) 
and the third-order essentially nonoscillatory (ENO3), the speed of the MP 
algorithm is of high rate, which makes the MP algorithm an efficient and 
effective procedure to solve the homogeneous conservation laws [5]. It is 
well-known that the MP scheme has higher accuracy and less dissipation 
compared with the mentioned schemes [1, 6, 5]. Moreover, the high-order 
WENO methods may not have monotonicity property but the MP scheme 
preserves monotonicity property [1]. Therefore, due to the mentioned reasons, 
we consider the MP method to discretize the homogeneous conservation laws 
part in the splitting process. 


4 Main results 


4.1 Monotonicity preserving 


In this section, we determine results for the basic splitting (2), (3) and the 
Strang splitting (4)—(6) with the above mentioned MP process for the homo- 
geneous conservation laws part. These schemes are called the MP basic split- 
ting scheme and the MP Strang splitting scheme, respectively. We state the 
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following theorem which determines the condition of monotonicity-preserving 
on the MP basic splitting scheme. 


Theorem 1. The MP basic splitting scheme is monotonicity-preserving un- 
der the CFL condition v < ths where v is CFL number. 

Proof. Suppose that the sub-problem (2) is computed with the TVD-RK3 
method in time and the MP scheme of five-order in space. Therefore, we 
suppose that uj—1 < fj_1 < uj, and uj < fj,1 <u". For increasing data, 
these two assumptions imply that (w7’ — uj—1) < (a +1)(uj — uj-1). After 


applying the relation (uY” — uj_1) < (a +1)(u; — uj—1), we will have 


0< (Fiza — fya) < (uw? — uj-a) S$ (@+ 1)(uy — uy-1). (8) 


1 1 
2 2 


Multiplying —v and then adding u; to the relation (8), we will obtain 


(l—v(l+a))u; +v(a+ ljuj-1 < u* < u;, 


where u* = uj; —v(fipa —fj-1) is the first step of the TVD-RK3 method. To 
hold uj-1 < u* < uj, the inequality uj-1 < (l—v(1+a))uj +u(at+ 1)uj-1 
must be satisfied, which in turn yields 


v(a + 1)(uj — uj—-1) < (ug — uy—1). 


Therefore, we have v(a +1) < 1. Since the second and third steps of the 
TVD-RK3 method are combinations of the first step and u; with positive 
weights, the final step u*(tn+1) is also monotonic. Consequently, the first 
step of the basic splitting scheme under the CFL condition v(a +1) < 1 is 
monotonicity-preserving. The sub-problem (3) is computed with the TVD- 
RK3 method with initial condition u**(t,) = u*(tn41) and also the TVD- 
RK3 method is monotonic; It is well known that if sub-problems is solved 
by a high resolution or a monotonicity-preserving method, then the split- 
ting methods (such as the basic splitting, Strang splitting and dimension 
splitting schemes) are monotonicity-preserving and convergent to the correct 
solution of original equation [4, 1, 29]. however, the MP basic splitting is 


monotonicity-preserving under the CFL condition v < ia: 


Here, we present the following theorem to estimate the error of the MP 
basic splitting method. 


Theorem 2. Suppose (2), (3) are computed with TVD-RK3 in time and the 
MP process of five-order in space. Then, in the absence of the splitting error, 
the global order of the error can be introduced as O(At?) + O(Az’). 


Proof. The computation of equation (2) with initial data un yields approx- 
imations u*(tn41) with an order error of the form O(At®) + O(Az®) [5]. 
More pointedly, the evaluation of equation (3) with TVD-RK3 method would 
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result in local error of O(At?) in time [36]. Consequently, Due to the 
stability of sub-schemes the global order of the error can be obtained as 
O(At?) + O(Ax°), when the splitting error is neglected. 


Also, to make sure that the convergency of the MP basic splitting method 
is held, we state the following lemma. 


Lemma 1. The MP basic splitting scheme is convergence to the true solution 
of original equation (1). 


Proof. As was mentioned above, subproblem (2) is computed with a high- 
resolution method in space (i.e., the MP process of five-order) and the method 
of TVD-RK3 in time with initial condition u*(t,) = un in interval [tn, tr41], 
while the subproblem is solved in interval [tn,t,41] only by the TVD-RK3 
method with monotone initial data u**(tn) = u*(tn41). Therefore, because 
sub-problems are solved using high-resolution methods (we know that high- 
resolution methods converge to the correct solution of the problem [29]) and 
the use of the well-known basic splitting method, also, because there is consis- 
tency between the approximation of the splitting method (2), (3) with the so- 
lution of original equation (1), the convergence of basic MP splitting method 
to the true solution of original equation (1) is guaranteed [4, 1, 29, 4]. 


4.2 Positivity preserving 


Next, we will investigate positivity property for the basic splitting (2), (3) 
sophistically in detail. We consider the case that the solutions u(x,t) of the 
basic splitting (2), (3) are nonnegative for all times t > 0. 

In the following, we state a lemma for positivity property of the basic 
splitting (2), (3) obtained using the MP method and the TVD-RK3 method. 


Lemma 2. The MP basic splitting scheme will be positivity-preserving under 
the CFL condition vy < Ta: 


Proof. Suppose that for t > 0 the problem (2) is solved with initial condition 
u*(tn) = Un > 0 in interval [t,,tn41] using the MP process in space and 
the TVD-RK3 method in time. Therefore, the solution of the MP method 
is monotonicity-preserving under the CFL condition v < its: However, the 
problem (3) is solved in interval [tn,tn+4i] by the TVD-RK3 method with 
monotone initial data u** (tn) = u*(tn41). Since the TVD-RK3 method with 
positive weights is positivity-preserving, therefore, u”*! := u**(tn11) > 0. 
As a result, the MP basic splitting is positivity-perserving under the CFL 


ees 1 
condition v < ro 


It is interesting to note that, the above-mentioned results are also valid 
for the MP Strang splitting method, i.e., the MP Strang splitting method 
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is convergent to the true solution of original equation (1) and also preserves 


monotonicity and positivity under the CFL condition v < a 


4.3 The case for diffusion source term 


In this section, we study the MP splitting process for balance laws with 
diffusion source term in the following form 


Ut + f(u)x = 1) Ura, (9) 


where f(u) is the flux function, and 7 denoting the diffusion coefficient. By 
applying the basic splitting method, the equation (9) is divided into two 
sub-problems as follows: 


uz; + f(u*)e = 0, for ty <t<tn4i, u*(z,tn) = u(z,tn), 


(10) 
(11) 


where u* and u** stand for concentration in the case of homogeneous conser- 
vation law and the diffusion process, respectively. Executing the MP process, 
equation (10) can simply be solved by using the initial condition of equation 
(9). The solution of this problem can be regarded as an initial condition for 
equation (11). Using the Crank-Nicolson method [8, 23], equation (11) can 
be numerically solved through an initial condition derived from the previ- 
ous step. Consequently, the solution to the original problem (9) is obtained. 
It turns out that the proposed process can stably and consistently estimate 
the true solution of the problem (9). This process can, however, be ap- 
plied to solving a wide variety of balance laws with possibly nonlinear source 
terms. As it is well-known, the effective and applicable strategy for solving 
the considered problem turns out to be implemented by the MP process. 
This method is a robust and powerful technique for treating balance laws. 
Moreover, it may be used in investigating problems involving diffusion source 
terms. To shed light on this issue, we will present the following remark which 
is, in essence, similar to the one given by Lemma 1. 


2k 2K 


up =nue., for th <t<tn4i, w*(a,tn) = u* (2, tn), 


Remark 1. It is notable that the proposed method is not only able to 
treat problems in the presence of reaction source terms, but also robustly 
and rigorously reconstructs the true solution as some diffusion source terms 
are imposed. Indeed, the proposed method can therefore retrieve the true 
solution of the original equation (9) by means of splitting method like what 
is mentioned in Lemma 1. This consideration makes us assert that this 
method can actually guarantee the convergency of the solution to the original 
equation (9). 
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5 Numerical results 


In this section, some test problems are provided to show the efficiency and 
effectiveness of the MP splitting schemes. In all tests, a = 4 and pp = 107 1° 
are considered. Meanwhile, all tests were performed by MATLAB routine 
codes. 


Test 


Test 


1. Consider the following equation 


uz + (au), = —Bu, a=land@=1,0<2<1, 
with initial condition 
u(x) = sin* (12), 


and the periodic boundary conditions, u(0,t) = u(1,t), at final time 
ty = 1. The results of the MP basic splitting scheme and the MP 
Strang splitting scheme are demonstrated with ale 0.05 and at =0.4 
in tables 2 and 3, respectively. The numerical results confirm the good 
efficiency and accuracy of the MP splitting schemes for this test. It 
is worth emphasizing that the numerical results for large CFL number 


a = 0.4 are also satisfactory. 
2. We solve 


ty + (au), = —Bu, a=land6=10, 0<t< 


Nl rR 


with initial condition given by 
uo(z) = 1000 cos (rx)’, 


and the periodicity conditions in the boundaries on the interval [0, 1], 
at final time tr = 5 and discrete points N = 40. Numerical solution of 
the advection-reaction equation with the MP basic splitting method for 
ae = 0.4 is depicted in Figure 2. From this figure, one can deduce that 
the proposed method is successful in retrieving the exact solution in a 
careful manner. Moreover, in comparison with the third order TVD 
basic splitting scheme [1], our proposed method is not only able to re- 
cover the true solution of the problem accurately, but also it preserves 
the accuracy of the numerical solution at extremum points. This situ- 
ation also happens for other test problems mainly because the spatial 
part of the homogeneous problem is discretized by the MP-process. Yet 
the lack of such property in the third order TVD-scheme makes its effi- 
ciency unsatisfactory. Furthermore, we consider the advection-reaction 
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equation as follows: 
+ (ou). =—Bu, a=land fs =1, (12) 


with initial condition 


uo(a) = (sin (x2)) cap (13) 


and the periodicity conditions in the interval [0, 1], at final time tr = 
2 and discrete points N = 64. Numerical solution of the advection- 
reaction equation (12)-(13) with the MP basic splitting method for 
cer = 0.2 is illustrated in Figure 1. In this test, we can deduce that the 
numerical solution of the MP basic splitting method, which preserves 
nonlinear stability such as monotonicity and positivity. To be more 
precise, this test is given specifically to confirm the positivity property 
of the MP basic splitting for CFL condition v < tis (as shown in 
Lemma 2) and it shows the MP basic splitting method is monotonicity- 


preserving for this CFL condition (as shown in Theorem 1). 


0.14 T T T T T T T T T 


—©— Numerical solution 
0.12 F Exact solution 4 


u(x,t) 


Figure 1: Numerical solution of the advection-reaction equation (12)-(13) with the MP 
basic splitting scheme at final time tf = 2 and N = 64, at = 0.2. 


Test 3. Consider the following equation governed by 
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ut + (au) =-—(6u, a=landé=1, 
with the initial condition 
ug(x) = sign(x), 


and the periodic boundary conditions in the interval [0,2]. The results 
are shown at final time ty = 1 in Figure 6. We get a good behavior for 
the MP Strang splitting scheme (even for CFL= 0.4). It is obvious that 
the MP Strang splitting scheme produces sharp results in the vicinity 
of the shock. 


Test 4. Let us take the equation 
ut + (au) = u(1 —u), a=1, 
with the initial condition 
ug(x) = exp (— 0.0127), 


and the periodic boundary conditions in the interval [0,2]. The results 
of the MP Strang splitting scheme are illustrated in Figure 4 with Bl a 
0.4, at final time tf = 1 and N = 100. However, for the final time 
tr = 1 and discrete points N = 100, the MP Strang splitting scheme 
is compared to the nonstandard finite difference (NSFD) scheme [33], 
which is stable and also is of first-order in time and space as follows: 


1 
n+l _ un ye 


J De ae 28 jot _ yn tiga), 


g(At) °  g(Az) a 


where $(h) = exp(h) — 1, and At = Az is an imperative condition to 
hold this scheme. We find that the MP Strang splitting scheme has a 
good resolution and convergence. Moreover, the result obtained by the 
MP Strang splitting scheme is free of the spurious oscillations. 


U 


Test 5. We apply the MP Strang splitting scheme to the equation 


ue + (5). =u(l-u), (14) 


with the initial value 


0.1+0.1sin (Qian); if O0<a<1l, 


0.1, else, 
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and the periodic boundary conditions in the interval [0,2]. The nu- 
merical solution is depicted at tf = 1.5 in Figure 5. We can see that 
the MP Strang splitting scheme can solve high-accuracy and resolution. 
Obviously, the dissipation of the MP Strang splitting scheme is much 
small. 


—©— Numerical solution 
Exact solution 


Numerical solution 
Exact solution 


Figure 2: Numerical solution of the advection-reaction equation with the MP basic 
splitting scheme (left), the third order TVD basic splitting scheme (right) at final time 
tp = 4 and N = 40, 3* =0.4 in Test 2. 


Table 2: Estimated L-errors and Lg-orders for advection-reaction equation with the 
MP basic splitting scheme in Test 1. 


N At = 0.05 at =04 
Lo-error Ly-order Lyo-error Ly-order 
20 6.13e—04 s 7.64e—04 - 
40 2.06e—05 4.89 5.43e—05 3.96 
80 6.58e—07 4.96 4.98e—06 3.45 
160 2.15e—08 4.94 5.64e—07 3.14 


Test 6. In this test, we apply the MP basic splitting scheme mentioned in 
subsection 4.3 to the equation 


1 
Ut + (5u’) =1Urx, 1=1, (16) 


with the initial condition 


uo(x) = sin(a), (17) 
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Table 3: Estimated L2-errors and L2-orders for advection-reaction equation with the 
MP Strang splitting scheme in Test 1. 


N &t = 0.05 Ri = 0.4 
Lo-error Lo-order Lo-error Lo-order 
20 6.10e—04 - 6.42e—04 - 
40 2.04e—05 4.90 2.46e—05 4.70 
80 6.50e—07 4.98 1.18e—06 4.38 
160 2.04e—08 4.99 8.8le—08 3.80 


—©— Numerical solution 
Exact solution 


u(x,t) 


-0.6 1 1 1 1 1 1 1 1 1 


Figure 3: Numerical solution of the advection-reaction equation with the MP Strang 


splitting scheme at final time t, = 1 and N = 200, at = 0.4 in Test 3. 


and the boundary conditions u(—4,t) = u(4,t) = 0 in the interval 
[—4, 4]. The results of the MP basic splitting scheme are illustrated in 
Figure 6 with Al = 0.4, at final time ts = 4 and N = 100. We can 
deduce that numerical solution given by the MP basic splitting scheme 
nicely converges to the true solution. 


Test 7. Now we consider the one-dimensional modified Burger’s equation 
(MBE) in the following form 


up +uP(ule — Nee =0, O<e<1, t>1, B=2, (18) 
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with the initial condition 


u(x, to =1) = 5, t>1,0<@<1,0<c<1,_ (19) 


1+te™ 


and the boundary conditions 


u(0,t) = u(1,t) = 0, 
Uz (0, t) = Uz (1, t) = 0, 
Ure (0,t) = Ure (1,t) = 0 


To apply the MP basic splitting, we split the original equation (18) into 
two the sub-problems P; and P»2 as follows: 


yortt 


Pr: uet (zr), = 9, 


Po: Ut = NUcze- 


Then, we get the solution of the MP basic splitting scheme mentioned in 
subsection 4.3 in the interval [0, 1]. The results of the MP basic splitting 
scheme are illustrated in Table 4 with At = 0.01, Ax = 0.0625, 7 = 
0.01, c= 0.5, at final time ty = 2. From Table 4 with these parameters, 
we can see that the MP basic splitting scheme error norms L. and 
Lg are in general better than the proposed methods in references [34, 
40]. MATLAB code of the MP process for this test is placed in the 
Appendix. 


Table 4: A comparison of the MP basic splitting scheme error norms Loo and Lz with 
the proposed methods in references [34, 40]. 


Schemes Loo-norm Lo-norm 

Sat[40] 7.4741e — 04 3.2033e — 04 
7.5978e — 04 3.4748e — 04 

MP basic splitting scheme 2.2677e — 04 1.4944e — 04 


6 Conclusion 


In this paper, the MP splitting schemes for solving balance laws are proposed. 
The MP splitting schemes are introduced to preserve the monotonicity and 
positivity properties of the numerical scheme for balance laws with linear and 
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0.995 + 


u(x,t) 
u(x,t) 


0.99 F 


—©— Numerical solution 


——©— Numerical solution 
1 1 0.985 1 H i 
0 


0.985 1 1 1 1 
0 


Figure 4: Numerical solution of the advection-reaction equation with the MP Strang 
splitting method (left), the NSFD method (right ) at final time ty = 1 and N = 100 in 
Test 4. 


nonlinear source terms. In these methods, the homogeneous conservation 
laws step is discretized by the fifth order MP method, and the reaction step 
and time integration are discretized by the third-order TVD-RK3 method. 
The numerical results for balance laws as well as the case for diffusion source 
term are also provided results that verify the applicability and efficiency of 
the proposed schemes. It is worth studing the MP procedure to the interface 
problems that have a wide application in science and engineering [11]. 
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APPENDIX: MATLAB code of the MP process 


function [fplus] = MP(u) 


Bi = 1/60; 
B2 = 4/3; 
alpha =4.0; 


eps = 10°-10; 
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0.6 T T T T T T T T T 


—©— Numerical solution 
Reference solution 


Figure 5: Numerical solution of the problem (14)-(15) with the MP Strang splitting 
method at final time ty = 1.5 and N = 100, at = 0.4 in Test 5. 


fhat =B1*(2 *u(1)-13*u(2)+47*u(3) +27*u (4) -3*u(5)) ; 
uMP = u(3)+DMM((u(4)-u(3)) ,alpha*(u(3)-u(2))); 
if ((fhat-u(3))*(fhat-uMP) <=eps) 


fplus=fhat ; 
else 
DJM1 = (u(1)-2*u(2)+u(3)); 


DJ (u(2)-2*u(3)+u(4)) ; 

DJP1 = (u(3)-2*u(4)+u(5)) ; 

DM4JPH = DM4(4*DJ-DJP1,4*DJP1-DJ,DJ,DJP1); 
DM4JMH = DM4(4*DJ-DJM1,4*DJM1-DJ,DJ,DJM1) ; 
uUL = u(3)+alpha*(u(3)-u(2)) ; 


uAV = 0.5*(u(3)+u(4)); 
uMD = uAV-0.5*DM4JPH; 
uLC = u(3)+0.5*(u(3)-u(2) ) +B2*DM4JMH ; 


umin = max(min(min(u(3) ,u(4)),uMD) ,min(min(u(3) ,uUL) ,uLC)) ; 
umax = min(max(max(u(3) ,u(4)) ,uMD) ,max(max(u(3) ,uUL) ,uLC)) ; 
fplus = fhat+DMM(umin-fhat ,umax-fhat) ; 

end 

return 

function A=DM4(W,X,Y,Z) 
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u(x,t) 


-0.04 —©— Numerical solution 
Reference solution 
-0.05 : 1 1 1 
4 -3 -2 1 0 1 2 3 4 


Figure 6: Numerical solution of the problem (16)-(17) with the MP basic splitting 
method at final time tr = 4 and N = 100, at = 0.4 in Test 6. 


A=0.125*(sign(W)+sign(X))*abs((sign(W)+sign(Y) )*(sign(W)+sign(Z)))* 
min(min(abs(W) ,abs(X)) ,min(abs(Y) ,abs(Z))); 
return 
function B=DMM(X,Y) 
B = 0.5*(sign(X)+sign(Y))*min(abs(X) ,abs(Y)); 
return 
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